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Swarming is an essential part of honeybee behaviour, wherein thousands of bees cling onto each 
other to form a dense cluster that may be exposed to the environment for several days. This cluster 
has the ability to maintain its core temperature actively without a central controller, and raises the 
question of how this is achieved. We suggest that the swarm cluster is akin to an active porous 
structure whose functional requirement is to adjust to outside conditions by varying its porosity to 
control its core temperature. Using a continuum model that takes the form of a set of advection- 
diffusion equations for heat transfer in a mobile porous medium, we show that the equalization of 
an effective "behavioural pressure" , which propagates information about the ambient temperature 
through variations in density, leads to effective thermoregulation. Our model extends and generalizes 
previous models by focusing the question of mechanism on the form and role of the behavioural 
pressure, and allows us to explain the vertical asymmetry of the cluster (as a consequence of buoyancy 
driven flows), the ability of the cluster to overpack at low ambient temperatures without breaking 
up at high ambient temperatures, and the relative insensitivity to large variations in the ambient 
temperature. Finally, our theory makes testable hypotheses for how the cluster bee density should 
respond to externally imposed temperature inhomogeneities, and suggests strategies for biomimetic 
thermoregulation. 



I. INTRODUCTION 

Honeybees are masters of cooperative thermoregula- 
tion, and indeed need to be able to do so for surviv- 
ing during winter, raising their brood, frying predatory 
wasps, or during swarming. Swarming is an essential part 
of colony reproduction during which a fertilized queen 
leaves the colony with about 2,000-20,000 bees, which 
cling onto each other in a swarm cluster, typically hang- 
ing on a tree branch, for times up to several days, while 
scouts search for a new hive location [1 j. During this 
period, the swarm cluster regulates its temperature by 
forming a dense surface mantle that envelopes a more 
porous interior core. At low ambient temperatures, the 
cluster contracts and the mantle densities to conserve 
heat and maintain its internal temperature, while at high 
ambient temperatures the cluster expands and the man- 
tle becomes less dense to prevent overheating in the core. 
Over this period, when the swarm is in limbo before mov- 
ing to a new home the cluster adjusts its shape and size 
to maintain the bees at a tolerable temperature, and is 
able to regulate the core temperature to within a few 
degrees of a homeostatic set point of 35° C over a wide 
range of ambient conditions. 

The swarm cluster is able to control its core temper- 
ature without a centralized controller to coordinate be- 
haviour in the absence of any long-range communication 
between bees in different parts of the cluster [2 j. In- 
stead, the thermoregulation behaviour of a swarm clus- 
ter emerges from the collective behaviour of thousands 
of bees [3] who know only their local conditions. So how 
can bees in a cluster, each acting on very limited local 



information, control a core temperature that they are ig- 
norant of? Early work on swarm clusters, and the related 
problem of winter clusters [H [5] , used continuum models 
for variations in bee density, and temperature as deter- 
mined by the diffusion of heat in a metabolically active 
material. Most of these models [6H9] assumed that the 
bee know their location and the size of the cluster, con- 
trary to experimental evidence. However, a new class of 
models initiated by Mverscough [TP] are based on local in- 
formation, as experimentally observed; these models are 
qualitatively consistent with the presence of a core and 
mantle, but are unable to maintain a high core temper- 
ature at low ambient temperatures. Further refinements 
of these models that account for bee thermotaxis and 
also use only local information[TT| [T2] , yield the observed 
mantle-core formation and good thermoregulation prop- 
erties at low ambient temperatures, while allowing for an 
increased core temperature at very low ambient temper- 
atures, observed in some clusters [2} [T3HT5] . However, the 
thermotactic mechanism that defines these models of bee 
behaviour causes the cluster to break up at moderate to 
high ambient temperatures, unlike what is observed. 

Here we present a model for swarm cluster thermoreg- 
ulation that results from the collective behaviour of bees 
acting based on local information, yet propagates infor- 
mation about ambient temperature throughout the clus- 
ter. Our model yields good thermoregulation and is con- 
sistent with experiments at both high and low tempera- 
tures, with a cluster radius and mantle density qualita- 
tively similar to observations, without leading to cluster 
breakup. In Section 2, we outline the basic principles and 
assumptions behind our model. In Section 3, we formu- 
late our model mathematically, characterize the number 
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of parameters in it, and solve the governing equations us- 
ing a combination of analysis and numerical simulation 
in Section 4, and compare our results qualitatively with 
observations and experiments. We conclude with a dis- 
cussion in Section 5, where we suggest a few experimental 
tests of our theory. 



II. MODEL ASSUMPTIONS 

Any viable mechanism for cluster thermoregulation 
consistent with experimental observations should have 
the following features: a) The behaviour of a cluster must 
result from the collective behaviour of bees acting on lo- 
cal information, not through a centralized control mecha- 
nism, b) The bee density of a cluster must form a stable 
mantle-core profile, c) The cluster must expand (con- 
tract) at high (low) ambient temperatures to maintain 
the maximum interior "core" temperature robustly over 
a range of ambient temperatures. 

This suggests that any quantitative model of the be- 
haviour of swarm clusters requires knowledge of the 
transfer of heat through the cluster, the movement of 
bees within the cluster, and how these fields couple to 
each other. The basic assumptions and principles behind 
our model are as follows: 

1. The only two independent fields are the packing 
fraction of bees in the cluster in the cluster p 
(equivalent to density) and the air temperature T. 
These then determine the bee body temperature 
and metabolism which can be written as a function 
of the local air temperature, while convection of air 
in the form of upwards air currents depends entirely 
on the global bee packing fraction and temperature 
profiles. 

2. We treat the cluster as an active porous struc- 
ture, with a packing fraction-dependent conductiv- 
ity and permeability. Bees metabolically generate 
heat, which then diffuses away through conduction 
and is also drawn upwards through convection. The 
boundary of the cluster has an air temperature that 
is simply the ambient temperature. 

3. Cold bees prefer to huddle densely, while hot bees 
dislike being packed densely. In addition, bees at- 
tempt to push their way to higher temperatures. 
The movement of bees is determined by a be- 
havioural variable which we denote by "behavioural 
pressure" P5 (p, T), which we use to characterize 
their response to environmental variables such as 
local packing fraction and temperature. In terms 
of this variable, we assume that the bees move 
from high to low behavioural pressure, a notion 
that is similar in spirit to that of "social forces" 
used to model pedestrian movements [16]. Here, we 
must emphasize that behavioural pressure is not a 



physical pressure. For packing fraction to be un- 
changing, behavioural pressure must be constant 
throughout the cluster. 

4. We assume that the number of bees in the cluster 
is fixed and that the cluster is axisymmetric with 
spherical boundaries. At higher temperatures, the 
clusters become elongated and misshapen, so that 
the assumption is no longer accurate; however, this 
does not change our results for thermoregulation 
qualitatively. In general, to determine the shape 
of the cluster, we must account for both heat and 
force balance, but we leave this question aside in 
the current study. 

A model based on these assumptions can be used to 
study both the equilibria and dynamics. However, sys- 
tems which reach equilibrium quickly are best under- 
stood in terms of their equilibrium behaviour rather than 
the specifics of how this equilibrium is reached. Since a 
swarm cluster is a constantly changing network of at- 
tachments between bees, as bees grab onto and let go 
of nearby bees, and can also detach and reattach them- 
selves at different points on the surface of the cluster, 
these "microscopic" dynamical processes allow the clus- 
ter to quickly equilibrate to changes in ambient condi- 
tions [21 [15] • Thus, we will mostly focus on the resulting 
equilibria, but consider the slow dynamical modes that 
allow it to respond to large scale weak forcing, as they are 
particularly important in the context of cluster stability. 

III. MATHEMATICAL FORMULATION 

The two independent fields in our model are bee pack- 
ing fraction p(r, t) G [0,1], and the air temperature 
T(r, t)\ while both of these are functions of space and 
time, we will focus primarily on the equilibrium be- 
haviour of these fields. For a static cluster of bees mod- 
eled as an active porous medium that generates heat 
and is permeable to air, the heat generated metabolically 
must balance heat lost due to conduction and convection, 
so that 

P M(T) + V • (k(p)VT) - Cu ■ VT = | ?en 

T = T a \? e sn, 

where k(p) is the packing fraction-dependent conductiv- 
ity, M. (T) is the metabolic heat production rate of the 
bees per unit volume, and C is the volumetric heat ca- 
pacity of air. We model the conductivity of the cluster as 
arising from a superposition of random convection cur- 
rents within the cluster which are suppressed at high bee 
packing fraction and the bare conductivity of the bees 
treated as a solid, and approximate this by a function 
k(p) = ko^^. Although this form diverges as p — » 
and random convection currents are unsuppressed, the 
p never vanishes in the interior of the cluster, so that 
this limitation is not a problem. Likewise, by bounding 
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p from above, we prevent k from vanishing in the cluster. 
We further assume that the mean flux per unit area u is 
determined by Darcy's law for the flow of an incompress- 
ible buoyant fluid through a porous medium, so that: 

u = [gap air (T -T a )z- VP] • k( P )/t] | ?Gn (2) 
V • u = 0, P = y e6Q . (3) 




FIG. 1: A cluster at an ambient temperature of a) 11°C 
and b) 27° C, approximately 12 cm across (Photo 
courtesy of [17] ) . Note that number of bees inside the 
cluster is nearly constant; change in cluster width is a 
result of changes in bee packing fraction pQ. c) 
Schematic of interior (11), and boundary (£0), with 
mantle-core structure, s, z are radial, vertical directions 
in polar coordinates. 



Here Eq. Q relates the flow velocity to the effects 
of thermal buoyancy and the pressure gradient, while 
Eq. ([3| is just the incompressibility condition, with n (p) 
the packing fraction-dependent permeability, a is the co- 
efficient of thermal expansion, p a i r is the density of air, 
r] is the viscosity of air, g is acceleration due to gravity, 
and P is air pressure. We assume that the permeabil- 
ity of the cluster may be approximated via the Carman- 

Kozeny equation n(p) = k,q ^"f^ , used to describe the 
permeability of randomly packed spheres [T8] . 

In general, the bee metabolic activity is not a con- 
stant, and depends on a number of factors such as tem- 
perature, age, oxygen and carbon dioxide concentration 
etc.[2j [19j [20]. To keep our model as simple as pos- 
sible, we start by assuming that the metabolic rate is 
temperature-independent, with M(T) = Mo and show 
that this is sufficient to ensure robust thermoregulation, 
setting apart the details of the calculations that show 
that our model can also yield robust thermoregulation 
with a temperature-dependent metabolic rate (Append. 

To close our set of equations, we still must relate p(r) 
to T(r), which we do by making a hypothesis that bees re- 
spond to local packing fraction and temperature changes 
by changing their packing fraction to equalize an effective 
behavioural variable which we call the behavioural pres- 
sure. Our formulation of this behavioural pressure relies 



on two assumptions that are based on observations: 

1. In their clustered state, bees have a natural packing 
fraction which is a function of the local tempera- 
ture. This natural packing fraction decreases with 
increasing temperature, and increases with lower 
temperature, until it reaches a maximum packing 
fraction pmax- Effectively, cold bees prefer to be 
crowded, while hot bees dislike being crowded, con- 
sistent with a variety of experiments in the field and 
in the laboratory [3 US HZ| . 

2. In addition to having a temperature-dependent 
natural packing fraction, bees also like to push 
their way towards higher temperatures. This will 
cause areas of equal local temperature to pack more 
densely at low ambient temperatures, consistent 
with observations [2]. 

With these constraints in mind, a minimal model for 
bee behavioural pressure suggests the piecewise linear 
function 

p r T \ = f -X-T+ |p-p m (T)| p<Pmax 

bK ' \ oo p> p maxj y ' 

where p m (T) = mm{p max , po — p'T} is the natural 
packing fraction, with the constants p' and x having 
units of temperature -1 . We emphasize that our model 
assumes that individual bees have no independent home- 
ostatic set points for temperature or packing fraction; 
the behavioural pressure depends on a combination of 
p and T. Furthermore, we note that the bee packing 
fraction in a cluster at equilibrium can never be lower 
than the natural packing fraction; in our formulation, 
the behavioural pressure is always higher at lower pack- 
ing fractions and indeed is the cause of the basic aggre- 
gation behaviour that creates the cluster. Additionally, 
the absolute behavioural pressure is of no consequence; 
only gradients and relative values are important. Further 
evidence for a behavioural pressure comes from experi- 
ments [4 who observed a relation between the ambient 
temperature and the core density, even when the core 
temperature remained approximately constant. 

To complete the formulation of our problem, we need 
to specify boundary conditions for the temperature and 
packing fraction fields. As we have stated earlier, the 
surface bees will be at the ambient temperature; further- 
more, since they can freely expand and contract to mini- 
mize their behavioural pressure, we assume that they will 
be at their natural packing fraction p m (T a ). 

Comparing our model with earlier models such as the 
Myerscough model [10] and the Watmough-Camazine 
model [TT], we note that both fit into this general 
behavioural pressure framework( Append. [B| ). However, 
our work differs from these studies in that we synthesize 
and generalize the implicit relation between bee be- 
haviour and their environmental variables in terms of the 
behavioural pressure, choosing a form for behavioural 
pressure which combines elements from both models 
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and is consistent with experimental observations. In 
addition, we account for the role of fluid flow and con- 
vection of heat in the cluster, which breaks the vertical 
symmetry, and is potentially relevant in heat transport 
at high Rayleigh number Ra = 9apt £ rC • (P — T a ) nR. 

Dimensionless equations: To reduce the number 
of parameters in our model, we make our equations di- 
mensionless. Noting that the total number of bees in 
the cluster is constant allows us to define a dimension- 
less total bee volume V = JJJ pdv/Vo, where Vo is the 
total bee volume of a typical cluster, i.e. the average 
volume of a bee times the average number of bees in a 
cluster. Upon setting a characteristic length scale to be 

the radius of a sphere of volume Vo, Ro = \J we 

write fff pdv = (47r/3)V. Scaling the temperature so 
that typical ambient temperature of 15° C — » 0, and the 
goal temperature of 35° C — » 1, we use the dimensionless 
variables: 

T-15°C T a - 15°C 

^ ? , •> J- a ^ 



= min{p +T(p' + X) ~ X^a, Pmax} 
= min {po+T -Co+Ta'Ci, pmax} , 



(11) 



20°C 

-> • K A ^ • (^0 , k 



20°C 
k Q (2Q°C) 

Mo^l, p' p' ■ (20°C) , X ^X-(20°C) 
and write the dimensionless form of Eqs. ([!]) to Q as: 

pM + V • (k(p)VP) - u • VT = (5) 
u = [(T - T a ) £ - VP] • «(p), V • u = 0. (6) 

k(p)=k -^A K (p)=^.Q—A- (7) 
T = T a y e sQj P = |j/ g jq (8) 

= f |p-p m (T)| p<pmax /gx 

\ oc p > /) max , 

p m (T) = min {p max , p - p'P} , (10) 



for the packing fraction and temperature profiles of the 
cluster which depend on the dimensionless parameters 

Pmax, P0, P', Xj T a, k , « , TV [H]. 

Information transfer through equalization of 
behavioural pressure: On the outer boundary of the 
mantle, the air temperature is equal to ambient temper- 
ature, so that minimizing the behavioural pressure re- 
quires that the packing fraction at the mantle will be 
the natural packing fraction at the ambient temperature 
Pm{T a ). At equilibrium, (4) then implies that the be- 
havioural pressure in the mantle and thus throughout 
the cluster will be — T a • x- This means that throughout 
the cluster, we may write the local bee packing fraction 
as a function of both the local temperature and the am- 
bient temperature, i.e P&(p(P, P a ), P) = —T a -\ which we 
solve to find 

p(T,T a ) = min {p m (T) + X -(T Pa) i Pmax} 



where we have made substitutions cq = p' + x, c\ = — X- 
Intuitively, the Co term characterizes the sensitivity of 
core temperature to ambient temperature, but the clus- 
ter cannot fully adapt at low T a through this term alone; 
adaptation at lower ambient temperatures requires the 
c\ term, which is eventually responsible for overheat- 
ing in the core at very low T a . In Fig. [2] we graph 
the packing fraction p(P, T a ) obtained by tracing con- 
tours of equal P5 which allows us to write the equilib- 
rium local packing fraction everywhere in terms of the 
conditions at the boundary. Bees respond to their local 
conditions and move accordingly, and these variations 
in packing fraction propagate information about ambi- 
ent temperature throughout the entire cluster without 
long-range communication. Although we have mapped 
Pfr(p, P) — >• p(P, T a ) for just one choice of behavioural 
pressure, our approach will work for any equation of state 
Pfr(p, P) which uniquely defines a stable packing fraction, 
i.e. with dPb/dp > 0. 




FIG. 2: Graphical representation of Eqs. ^ to showing 
how local packing fraction is obtained through behavioural 
pressure (Online version in colour). At low T a , bees on the mantle 
will pack at p = pmax (upper circle), and behavioural pressure is 
high throughout the cluster leading to higher interior packing 
fraction(upper solid line) and overpacking. At high T a , bees on 

the mantle will pack more loosely (lower circle), and low 
behavioural pressure throughout the cluster leads to low interior 

packing fraction (lower solid line). The solid lines run along 
contour lines of equal P&(p, T), as behavioural pressure must be 
uniform through the cluster. Coefficients used are 
po = .85, p' = .75, p m ax = -8, % = -3. 



The set of equations^ j6]j7j[8j and 11 ) with the bound- 
ary condition that the surface temperature of the cluster 
is the ambient temperature completes the formulation of 
the problem to determine p(r), P(?). 



IV. SIMULATIONS AND RESULTS 

While the boundary of the cluster is assumed to have 
spherical symmetry, the temperature and packing frac- 
tion fields inside do not have to have spherical symmetry 
owing to the effects of convection. However, the fields still 
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have cylindrical symmetry, and therefore we can repre- 
sent p(r),T(r) as p(s, z),T(s, z), where s is the distance 
from the central axis and z is the height. With this coor- 
dinate representation, we solve the governing equations 
§§|7|[8)and[ll]) in a spherically bounded domain using 
a simple discretization scheme with 30 values of s and 60 
values of z (Append. |C 2| ). 

Our choice of the dimensionless parameter ko = 0.2 
is constrained by experiments [19], while we estimate 
Ko = 1 though a simple calculation assuming the bees 
to be randomly packed spheres, and p m ax — -8? slightly 
higher than the maximum packing fraction of spheres, as 
bees are more flexible. (Append. A 2). However, the pa- 
rameters defining bee movement and behaviour, namely 
co,ei, and po are experimentally unknown. Guided by 
the general observation that physiological performance is 
often improved by changing parameters while basic mech- 
anisms remain unchanged, we optimize these parameters 
to achieve robust thermoregulation, i.e. the core tem- 
perature remains close to 1(35°C) over a range of scaled 
ambient temperatures T a G [—0.7, 0.8] corresponding to a 
real ambient temperature T a G [0,30]°C. For the choice 
of parameters po = .85, Co = .45, ci = .3, we find that 
within this wide range of T a and a factor of three in V, the 
dimensionless core temperature stays in the range .7 — 1.3 
corresponding to a real temperature range of 29— 41°C, 
while the core temperature itself increases monotonically 
with po in an analytically solvable way (Append. A3). 



Hot T c 



Cold T a 



Our simulations also capture the qualitative mantle- 
core structure of the cluster (Fig. [3| with a dense mantle 
surrounding a sparse core. We also find that at high am- 
bient temperatures, the cluster expands and the mantle 
thins, and at low ambient temperatures the cluster con- 
tracts and the mantle thickens. Furthermore, we find 
that the core temperature, which is also the maximum 
temperature of the cluster, is higher at low ambient tem- 
peratures resulting from "overpacking" , consistent with 
experiments of [2j [T3HT5] , and predicted earlier [TT . Fi- 
nally, we note that the temperature profile is vertically 
asymmetric due to convection, causing the point of max- 
imum temperature to rise above the geometric center of 
the cluster, as observed in experiments [2j [15] . 



V. DISCUSSION 

We have shown that it is possible to provide a self- 
organized thermoregulation strategy in bee clusters over 
a range of observed ambient temperatures in terms of 
a few behavioural parameters. Our theory fits into 
a broader framework for understanding collective be- 
haviour where the organism responds to the environment, 
but in doing so, changes it, and its behaviour until it 
reaches a steady state. Here, our model takes the form 
of a two-way coupling between bee behaviour and local 
temperature and packing fraction, quantified in terms 
of an effective behavioural pressure whose equalization 
suffices to regulate the core temperature of the cluster 




FIG. 3: Adaptation with temperature-independent 
metabolism(Online version in colour), a) Comparison of 
temperature and packing fraction profiles at high (.8) and 
low( — .7) ambient temperature, where the dimensionless total bee 
volume V is 1. b) Core temperature as a function of ambient 
temperature and total bee volume shows that as V changes the 
core temperature increases but adaptation persists over a range of 
T a (also plotted to guide the eye.) c) Cluster radius as a function 

of ambient temperature and total bee volume showing how 
clusters swell with temperature, consistent with experiment. For 

bee packing fraction, we choose coefficients of 
po = .85, c = .45, c\ — .3, Omax — -8. For heat transfer, we choose 
coefficients of ko = .2, k,q = 1. 



robustly. Although our choice of the form of the be- 
havioural pressure is likely too simplistic, it is consis- 
tent qualitatively with experimental observations, and we 
think it provides the correct framework within which we 
can start to quantify collective behaviour. Furthermore, 
our strategy might be useful in biomimetic settings. Our 
formulation for behavioural pressure shows that in a clus- 
ter, bee packing fraction should depend only on local 
temperature and the temperature of bees at the bound- 
ary, which effectively control the surface packing fraction. 
This observation provides an immediately testable pre- 
diction: a cluster may be "tricked" into overheating its 
core by warming the bees just below the surface, while 
exposing the surface bees to a low temperature to in- 
crease behavioural pressure. As pointed out earlier, ex- 
periments and observations of bee core temperature and 
bee packing fraction [4] are consistent with these ideas, 
although a direct experiment of this type does not seem 
to have been carried out. 

Our model adjusts well in response to changes in am- 
bient temperatures, but it doesn't have the same level 
of tolerance to different total bee volume that honeybees 
exhibit. We have used a continuum model where the bees 
on the surface exposed to the ambient temperature equi- 
librate to that temperature, but in reality the first layer 
of bees is hotter than ambient temperature and supports 
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a large temperature gradient driven by heat from inte- 
rior bees. This means that the surface bees feel an aver- 
age temperature higher than ambient due to the interior 
bees, and and we predict this should give a large cluster 
a lower behavioural pressure than a small cluster at the 
same ambient temperature. This should reduce the sen- 
sitivity of core temperature to total bee volume, and we 
have confirmed this in simulations (Append. |E| . 

We close with a description of some possible exten- 
sions of our study. Currently, our model ignores changes 
in shape of the cluster associated with force balance via 
the role of gravity, and the associated effects on ther- 
moregulation. In reality, the cluster is a network of con- 
nections between bees which changes shape and size due 
to a combination of mechanical forces and heat balance, 
and a complete theory must couple these two effects as 
well. 

Our heat balance equation and estimates (Append. 



A 2 ) suggest that while convective terms are responsible 



for the asymmetry in the temperature profile, they do not 
play an important global role in thermoregulation. How- 
ever, our calculation assumes a uniform packing which 
does not accurately represent the microscopic structure 
of the cluster; at high ambient temperatures, swarm clus- 
ters are observed to "channelize", where channels open 
up to ventilate. Studying a simple "behavioural pressure- 
taxis" dynamical law(Append. |f| , we find no linear in- 
stability that leads to channelization without an "anemo- 
taxis" mechanism, where behavioural pressure increases 
with |u|; instead we see only only one kind of linear in- 
stability which comes from the mathematical necessity of 
fixing cluster radius. This raises the question of whether 
there are anemotaxis mechanisms. If not, how can chan- 
nelization result from bee-level dynamics and mechanics? 
We have also neglected active cooling, which includes fan- 



ning, evaporative cooling( which can give up to ~ 50° C of 
cooling), diffusion of heat through diffusion of bees, and 
effects of oxygen and C0 2 [H 123 El ES] . Finally, we 
have also neglect any implications of bee age distribution, 
despite knowledge of the fact that younger bees tend to 
prefer the core and produce less heat, while older bees 
prefer the mantle and can produce more heat [2j [T4t [T7] . 
Accounting for these additional effects will allow us to 
better characterize the ecological and possibly even evo- 
lutionary aspects of thermoregulation. 

Thermoregulation is a necessity for a wide variety of 
organisms. When achieved collectively, individuals ex- 
pend effort at a cost that accrues a collective benefit. 
The extreme relatedness of worker bees in a cluster and 
near-inability to reproduce implies that the difference be- 
tween the individual and the collective is nearly nonex- 
istent, so that cost and benefit are equally shared. How- 
ever, many other organisms are faced with the "huddler's 
dilemma" [24 ; expending individual metabolic effort is 
costly, and benefits a group that is only partially related. 
Because genetic relatedness, metabolic costs, individual 
temperature and spatial positions are all easily measur- 
able, a collective thermoregulatory system is an ideal 
context in which to study the tangible evolution of co- 
operation and competition by building on our framework 
theoretically to account for competition and cooperation. 
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Appendix A: Units, parameter estimation, and dimensional analysis 



1. Table of quantities and associated units 



Quantity 


Symbol 


Units 


Typical Values 


Bee Packing Fraction 


P 


Dimensionless 


- 2. - .8 


"I— ^ "TV /T i TJ TJ • T~~\ , 

Bee Metabolic Rate 


M 


Power /Volume 


- .001 - .OlW/cnr 5 PSJH21 


Heat Conductivity 


k 


Power/ (Distance • Temperature) 


- .0004- .006W/(cm°C)? (Append. A 2) 


Permeability 


K 


Distance 2 


(.05cm) 2 ? (Append. 


A2 




Darcy Velocity 


U 


Distance/Time 


~ lcm/s? (Append. 






Air Heat Capacity 


C 


Energy /Volume -Temperature 


1.2 • 10- d J/(mL°C) 


Gravity 


9 


Distance/time 2 


981cm/s 2 


Air Density 


Pair 


Mass /Volume 


1.2 • 10" 3 gram/mL 


Coefficient of thermal expansion 


a 


Temperature -1 


1/300°C 


Air Viscosity 


V 


Pressure • time 


1.8 • 10 _4 gram/cm • s 


Behavioural Pressure 


Pb 


? 


? 


Thermotactic Coefficient 


X 


Behavioural Pressure/Temperature 


? 



Many sets of dynamical equations for bee movement will result in the same equilibrium where behavioural pressure, 
whose units depend on the set of dynamical equations used, remains constant. For example, a simple taxis model is 
one where: 



dp 
~dt 



V J, J 



would mean that behavioural pressure has units of Distance 2 /Time. A more complicated evaporation/condensation 
model would have the form 



dp(r) 
dt 



\r~r'\ 2 



. x e 



(A(f)-A(r)) 



where a is the evaporation and condensation radius and T is some transfer coefficient would give unites of 1/Time. A 
yet more complicated model involving mechanical compressibility or viscosity would have yet another set of dimensions 
for behavioural pressure. All of these models would, however, yield the same static solution. 



2. Estimation of heat conductivity, metabolic rate, and permeability 



For the metabolic rate and heat conductivity we use estimates from the experiments of Southwick [19 , where a 
cluster of 4250 bees(608 grams) is put into a set of roughly planar, parallel honeycombs, and the temperature profile 
and oxygen consumption are measured. The bees are roughly uniformly distributed with a bee packing fraction p of 
about .5, for which the metabolic rate is roughly uniform, and the temperature is well approximated with a parabolic 
profile, with a temperature of 34°C at the core, and 11°C at the edge 9.5cm from the core, parallel to the combs. The 
combs are very well insulating, so heat transfer occurs primarily in the two directions parallel to the combs. Then, 
within the cluster: 

x 2 + y 2 
9.5cm 2 ' 
_ 1.0°C 
cm 2 

The oxygen consumption rate is measured to be 6.5mL/min which gives a volumetric metabolism of .0035 • W/ (cm 3 ), 
assuming a bee specific weight of 1 gram/(cm 3 ), and an oxygen to energy conversion of 3.5 ^ g L ^ = -0012^^. 
This metabolic rate agrees well with the experiments of Heinrich [2]. We now have all the pieces to calculate the 



T « 34°C - 23°C • 



V 2 T : 



23° C 
(9.5cm 2 ) 



9 



conductivity using the conduction heat balance: 

k y 2 T + ^=o^k= i - 7 - io : 3w 

cm°C 

At p = .8, the maximum packing fraction we allow, the conductivity becomes close to the value of 2.4- 10 _3 W/ (cm°C) 
for fur and feathers [25] as suggested by Southwick which gives us some level of confidence in the functional form 
k°^7T we nave chosen for conductivity. 

To estimate ^o, we use the Carman-Kozeny equation [18 j. The average bee weighs about .14 grams, which cor- 
responds to a sphere of diameter ~ .65cm. In the absence of any detailed information about the bee structure in 
the cluster, we treat the cluster as a system of randomly packed spheres. It would be interesting to measure and 
better understand convective gas and heat transfer within swarm clusters. Using this diameter in the Carman-Kozeny 
equation, we find: 

D 2 

^° = 180 ~ (- 05cm ) 2 • 

A typical c luste r has about 10,000 bees, which which is about 1.4 kg, giving R = 7cm. Plugging these values 
in(Append. |A3[ ), we find dimensionless conductivity and permeability to be: 

k « .2, ft ~ 1. 



3. Dimensional analysis 



Our conditions for heat balance imply that: 



pM(T) + V • (k(p)VT) - Cu • VT = | ?en 

T = T a y e sn 

u = [gap air (T -T a )z- VP] • k(p)/t) y en 
V • u = 0, P = y eS n. 



while our equation for behavioural pressure reads: 



-X'T+\p- Pm(T)\ p < p max 

OO p> Pmax, 

pm(T) min {pmax, Po - p'T} 

Setting our unit of length to be Rq, we make the transformation V — >• V/Rq- Then our heat balance equations read: 



1 „ „ , 1 

Cu • 

T = T a \? e sn 



P M(T) + -5 V • (k(p)VT) - — Cu • VT = | ?€n 

±1 ±1q 



gap air (T-T a )z-^-VP 
±io 



v • u = o, p = o y eS n. 



On making the substitutions: 



T -> T - 15 ° C ^ T„- T " 15 ° C 



35°C-15°C' 35°C-15°C 
leads to the goal temperature of 35° C yielding a dimensionless temperature of unity, and a typical ambient temperature 



10 



of 15° C corresponding to a dimensionless temperature that vanishes. Then our system of equations becomes: 
P M(T) + 35 ° C - 2 15 ° C v • ( k (p)vD - 35 ° C ~ 15 ° C • Cu • VT = | ?en 



Rn 



U 



(35°C - 15°C) ag Pair (T -T a )z- — VP 

Ho 

Vu = 0, 



T = T a \? e sn 

P = \reSQ. 



P b (p,T) 



-x(35°C-15°C)T+|p-p m (T)| p</w 
p m (T) = min {p max , [po - p' • 15°C] - p'T[S5°C - 15°C]} , 



We divide all terms in the heat equation by the base metabolism Mo to yield: 

pM(T)/M + ^"T^ V • (k(p)VT) - 35 °^" D 15 ° C • Cu • VT = |* n 



M R 



u 



(35°C - 15°C) agp air (T -T a )z- —VP 

±1q 



T = T a y e sn 

V • u = 0, P = | ?G5n . . 



Making the substitution u — >> u 



35°C-15°C 
A4 Ro 



• c 



and the substitution 



P-^ 



(35°C- 15 o C)a 5/ 9 air R ' 



leads to the set of equations: 



pM(T)/M + KA „ 2 V • (k(p)VT) - u • VT = | ?etl 



35 °C-15° C. 

T = T a Ipg^n 

r/T tm- V7P1 (35°C-15°C) 2 o^w-C 
u = (T - T a z — VP /s(p) ?en 

V • u = 0, P = | Fe5n . . 

Finally, making the substitutions for the coefficients 



M R 2 

(35°C-15 C) 2 aff^ r -C 



Po^Po- p' (15°C) , p'^p'- (35°C - 15°C) , X -»• X • (35°C - 15°C) . 
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leads to our full dimensionless set of equations for heat balance: 



pM(T) + V • (k(p)VT) - u • VT = | FeC 

T = T a y e sn 

u=[(T-T a )z-VP]-K(p) | ?en 
V-u = 0, P = | ?€m . 



while the dimensionless behavioural pressure reads: 

n(p,T) = 



-X + IP - PmCOl P < Pmas 
OO p > p m axj 

pm(T) = min {p max , p - p'T} , 



Our model has seven parameters pmaxi Po, p' •> X? T a , ko, fto- 



a. ./Vote on further dimensional analysis 

We note that, when the metabolic rate is temperature independent, the goal temperature and the typical indepen- 
dent ambient temperature have no bearing on the actual behaviour of he model, only on whether it represents effective 
thermoregulation. Then, we may set the temperature where the packing becomes maximally dense, T pac k e d = pQ ~ p ^ ax ? 
to be zero, and the temperature at which the packing fraction becomes zero, T empty = ^7, to be 1. We can then make 
slightly different substitutions for the coefficients : 

M^l, k -> k • T ^Pty -^Packed 

(^empty ^packed) ^9 Pair ' C 



PO ^ Pmaxi P ^ Pmaxi X ^ X ' (^empty ^packed) i 

and find that there are now only five free parameters pmaxi X? T a , kg, Ko- We do not carry out this extra analysis in 
the main body of the paper because this causes us to lose sight of what the goal core temperature and typical ambient 
temperatures are. 



Appendix B: Behavioural pressure formalism and its antecedents in previous models 

To understand how our behavioral pressure formalism fits in with previous models, we compare them within this 
framework. The Myerscough model assumes that the bee density depends only on local temperature, and thus can 
be written as Pb(p, T) = \p — p m |, where p m = 8bees/cm 3 • (l — ^q). The Watmough- Camazine model defines a 
dynamical law for the bess density via the equations: 

p = -V • J 

J = -/x(p)Vp-px(T)VT, 

where fi(p) > is a motility function, and x(T) is a thermotactic function. This may be written as 

>(/>), 



-p- 



L P 



-Vp + X (T)VT 



-pVPb 
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where the behavioural pressure P5 is defined as: 

JMP,T)= j^dp' + jx{T')dT'. 

V v ' V v ' 

Density Component Temperature Component 

We note that as p{p) > for all p, the density component is minimized as p — » 0, and the density at the surface must 
be fixed to prevent the cluster from falling apart, unlike in our behavioral pressure framework. Additionally, because 
the behavioural pressure can be divided into temperature and density components, the point of highest density will 
always be at the same temperature, regardless of size or ambient temperature, which is not observed experimentally. 
This is in contrast with our formalism. 



Appendix C: Numerical methods 

1. Method of solving for equilibrium 

To reach equilibrium, we use an iterative scheme, described by the following pseudocode: 
T(r) <- T a 

max 

n (?\ ^ J P™x ■ |P| < R 

p[r) *~ \ : |r| > R, 
repeat 

(?) <- p(T(r) 
p(r) <- p(r) + c p (p new (r) - p(r)) 

Expand or contract the bee packing fraction and temperature profiles to normalize total bee volume 
7Z <— ? -jrr- Scaling ratio 



fffp 
R^R >K 

p(f) <- p(r/K) 
T(r) ±- T(f/K) 

Solve for u at fixed temperature and density, then solve for T at fixed M and u 
Find u such that: u = [(T — T a ) z - VP] • n(p), V • u = 
Find T new (r) such that: pM(T) = -V • (k(p)VT new ) - u • VT new 
T(r)^T(r) + (T new (r)-T(r)) 
until converged 

The intermediate steps can be solved as a system of linear equations. The solution is considered to be converged 
when p new = p, T new = T, 1Z = 1 to within 10 -10 . 



2. Discretization of space 



To solve for the temperature and density profiles, we must first discretize the system. While spherical symmetry 
is broken due to convection, the system retains rotational symmetry about its axis. We therefore use cylindrical 
coordinates, where each cell is given indices (i, j), and has coordinates which represent the distance from the central 
axis Sij and the vertical coordinate 2^, where 

n is the radius of the cluster in cells. All cells with s?- + z% < R are in the interior of the cluster, while all cells with 
s ij + z ij — ^ are a ^ ^ ne exterior of the cluster, subject the the boundary conditions — T a , = 0, pij = 0. 

The volume of each cell with coordinates is Vij = 2iTW 2 Sij, where w = ^ is the width of each cell. Each cell 
(i, j) neighbors four other cells, (i,j + 1), — 1), (i + 1, j), (i — 1, j), with the exception of cells which border the axis 
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(i = 0), which only have three neighbors. The area shared by cell and its outside horizontal neighbor (i + 1, j) 
is 2ttw (sij + The area shared by cell and its vertical neighbor (i, j + 1) is 2irwsij. 




FIG. 4: Cell layout of a cluster 30 cells in radius. Interior cells are coloured white, exterior cells are coloured grey. 



a. Heat equations 



For heat balance, we discretize our heat equations as 
dQij MiTi/lpijVij - ~~T~ ' H ( k (fti), Hpi>j>)) {Ti-r ~ To) + n >->r \-° Kw) 2ij + * ("".,,.'/) = °- 

Z Z 



Metabolism 



conduction 



convection 



The first term corresponds to metabolic heat generation, and the second term is heat conduction, where Ylu'j') ^ s 
the sum of all (i',j f ) neighboring with the heat conductance between two neighboring cells depending on the 

. The third term represents convective heat transfer, 



l/a+l/6 



harmonic mean of the conductance of each cell; H (a, b) 
with Uijj'j' the air flow from cell to (i f ,j'). 

We have chosen an "upwinding" scheme [26]; when air flows out of a cell, the outwards heat flux is determined by 
the temperature of that cell. When air flows into a cell from a neighboring cell, the inwards heat flux is determined 
by the temperature of the neighboring cell where the air originates. This scheme uses the Heaviside step function: 



£>(. 



x < 
x > 0. 



b. Solving for buoyancy driven flow 



The air flow from cell ij to neighboring cell i' j' is given by: 



A 



W 



(Pij Pi'j') H~ ( z i' 



Zij) 



T 



■ Tifjt 



where the air conductance between two cells again depends on the harmonic mean of the permeability of each cell, 
and the pressure is set so that air flow is conserved in every cell i.e. u ij,i'j' ~ f° r a ^ ce ^ s (hj)- This yields 

a set of linear equations that can be solved easily. 
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Appendix D: Temperature dependent metabolic rate 

To formulate our temperature-dependent metabolic rate, we note that bees have a higher base metabolic rate at high 
temperatures than at low temperatures. At moderate temperatures, bees on the mantle keep their body temperatures 
approximately 3°C above ambient temperature, and at below 15° C, they will "shiver" to keep their body temperature 
at 18°C[15]. Assuming a constant coefficient of thermal transfer between a bee and the surrounding air, this leads to 
a formulation of metabolism by shivering at air temperatures below 15° C, and gives us the full piecewise function for 
metabolic rate. Our metabolic rate for high temperatures comes from experiments involving oxygen consumption in 
swarm clusters [T5] . 

f 1 _L 15°C-T . rp iron 

where the base metabolism, Mo has units of power/ volume. 



6 
4 

M(T) 

2 




-1 -0.5 0.5 1 

T 

FIG. 5: Metabolic rate as a function of temperature. 

The simulated cluster with temperature-dependent metabolic rate has the same qualitative features as when we 
use temperature-independent metabolic rates. We set Mo to be half as high as when M is temperature-independent, 
because now it represents the minimal metabolic rate, not the uniform one (Fig. [6]). 



Hot T a T Cold T a Hot T a P Cold T a 




T a T a 



FIG. 6: Adaptation with temperature-dependent metabolic rate (Online version in colour), a) Comparison of 
temperature and packing fraction profiles at high (.8) and low (—.7) ambient temperature, where the dimensionless 
total bee volume V is 1. b) Adaptation is nearly the same as for temperature independent metabolic rate. T a is also 
plotted to guide the eye. c) Cluster radius is also nearly the same as for temperature independent metabolic rate. 
For packing fraction, we choose coefficients of po = .85, Co = .5,ci = .25,p ma;E = .8. For heat transfer, we choose 

coefficients of ko = .4, kq = 2. 
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Appendix E: Accounting for the role of finite bee size on thermoregulation 

In the paper, we have defined the temperature at the boundary to be the ambient temperature, and we assumed 
the surface bees feel the ambient temperature, which sets the behavioural pressure accordingly. In reality, bees point 
their heads inwards, and feel a temperature gradient driven by the heat produced by interior bees. Therefore, it may 
be more realistic that the behavioural pressure is set by the temperature a slight distance inwards from the surface. 
If we include the effects of convection, this implies that spherical symmetry must be broken, and the temperature 
becomes not just a function of distance, but also dependent on angle. Therefore, to close our set of equations without 
having to delve deeper into the question of cluster shape, we must neglect upwards convection and only consider 
conduction. This gives us the system of equations: 

pMo + V • (k(p)VT) = 0, k(p) = k • 

9 

where the behavioural pressure is now set by the temperature a distance of Lb ee inside the cluster: 

= min{/0o+T-Co+T(R-L b ee)-Ci, Pmax} • (El) 

Assuming Lb ee = 1 cm ? slightly shorter than the body length of a worker bee, for an average cluster radius of 7 cm 
we find that the dimensionless Lb ee of ~ .14. We simulate the system for dimensionless total bee volumes of .5, 1, 
and 3, with the same parameters of cq = .45, c\ = .3, pmax = -8 as were used in paper. In the first system, to test 
thermoregulation with this effect, we choose Lb ee = -14, and compensate for the slightly lower average behavioural 
pressure by increasing po to .95. In the second system, we test thermoregulation without this effect, and so we choose 
Lbee = 0,po = .85, as we did in the main body. We find that for large clusters, the temperature gradient created by 
the interior bees lowers behavioural pressure and loosens the cluster, mitigating overheating in the core and sensitivity 
to total bee volume. 



a) b) 




-0.5 0.5 -0.5 0.5 

Ta Ta 

FIG. 7: Comparison of core temperatures(Online version in colour), a) Behavioural pressure is set by the 
temperature slightly below the surface of the cluster, and sensitivity to total bee volume is mitigated, b) 
Behavioural pressure is set by ambient temperature, and overheating of the core is much worse. T a is also plotted in 

each to guide the eye. 



Appendix F: Results of linear stability analysis 

Solving for linear stability using simple "behavioural pressure-taxis" dynamics (Append. |F 1| ), we find that all 
clusters simulated at a temperature-independent metabolic rate are stable. However, at low ambient temperatures, 
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clusters with a temperature-dependent (Append. [D| metabolic rate can be linearly unstable via an overheating 
instability, which we believe to be a relic of fixing the boundaries of the cluster. In this instability, bees from the core 
move to the mantle, increasing the mantle thickness and insulation. This causes the core to heat up, increasing its 
behavioural pressure causing even more bees to move from the core to the mantle, leading to eventual runaway. We 
only see this instability when using a temperature-dependent metabolic rate, where an increased core temperature 
results in a greater net metabolic rate, aggravating the problem. 



Cluster 
Profile 




Instable 
T Eigenvector P 



w 



0.20.40.60.8 1 1.2 

s 



0.20.40.60.8 1 1.2 

s 



0.20.40.60.8 1 1.2 

s 




0.20.40.60.8 1 1.2 

s 



FIG. 8: a) Cluster profile, b) Contour plots of instable eigenvector of linear response matrix(Online version in 
colour). Circle represents boundaries of the cluster, and a temperature-dependent metabolic rate is is used. 
Jo = 1, V = 1.5, T a = —.7. po = .85, Co = .5,ci = .25,p ma;E = .8. For heat transfer, we choose coefficients of 
ko = .4, fto = 2. Note that the packing fraction very close to the boundary is fixed because p = p m ax- 



Dynamical behavior which allows the cluster radius to vary would suppress this instability, as bees moving from 
the core to the mantle would result in an expansion of the cluster, increasing the surface area and cooling the core. 
However, a dynamical model which allows the boundaries of the cluster to change requires a better understanding of 
the bee-level structure and the mechanics within a swarm cluster, and we leave this aside here. 



1. Methods for calculating linear stability 

Having solved for the equilibrium state, we want to find if this state is stable or unstable. To do so, we must 
define some dynamical laws for bee movement. Choosing a simple "behavioural pressure-taxis" behaviour, where 
J oc — pVPb, = — V • J gives us the complete set of dynamical equations: 

pf = pM(T) + V • (k(p)VT) - Cm • VT = | ?en 
u = [gap air (T -T a )z- VP] • K,(p)/rj 

J = -J • pVP b , J = \z e sn, 

Here we vary Jo over a wide range to reflect the large variations in bee movement and temperature time scales, and 
define F, G to be the functionals that determine the dynamics of the system. We also emphasize some issues about 
these choices: a) Our form assumes a substrate that the bees move on. In reality, the bees are the substrate, and 
changes on one side of the cluster propagate mechanically to other parts of the cluster at a rate faster than the taxis 
rate, b) Bee movement is not necessarily a local movement down pressure gradients. Bees can disconnect from the 
cluster and reattach at a different point, which doesn't fit into the local gradient picture, c) There is a discontinuity 
in behavioural pressure and packing fraction at the boundary of the cluster, so this taxis model does not explain 
how the boundary of the cluster can change. Within these limits then, this choice of behaviour gives us a full set 
of differential equations. To determine if the equilibrium state of the system is linearly stable, we must determine 



= pF\p,T\, T = T a y eSQ 
V • u = 0, P = | ?€m 

$ = -VJ y e n=G[p,T] 
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whether the linear response matrix: 



( dF dF\ 
% $6) 
dT dp J 

has positive eigenvalues. We note that the equilibrium temperature and bee packing fraction profile we solved for is 
symmetric with respect to rotation about the central axis((/> direction), and therefore we may partition the space of 
perturbations into subspaces defined by wave number k : AT(s, z, <j>) = e i ^5T(s, z), Ap(s,z,4>) = e i ^ <t> 5p{s, z). 
These temperature and bee packing fraction perturbations will in turn give a change in airflow, pressure, and be- 
havioural pressure: Au(s, z, <j>) = e ik ^<Ju(s, z), AP(s,z,</)) = e i ^ <t> 5P{s, z), AP h (s,z,<j>) = e^SP^s.z). All of 
these will change the temperature and bee packing fraction time derivatives which will be proportional to e* k< ^. For 
each wave number k^, we construct the stability matrix and study its spectrum. To first order: 



Tp 



-u • W sz + V sz • kV s 



k • | ^ 



AT+ 



[MAp + V sz k p V sz T]Ap - Au • V SZ T 



P = V sz • (pV sz AP b ) 



pAP b 



dP b A dP b . m 

where \7 SZ is the gradient in the s and z directions, Air — ^j 7 ^ , k p = • Regions where p = p ma x are locked 
at pmax an d not allowed to vary in bee packing fraction. 



a. Solving for Au, AP 

The above equations depend on Au, which is determined by: 

Au = k[ATz - (AP)] + Ap[(T -T a )z- VP]. 

Au must have the form: 



Au 



Au s , - ftV AP0 



Au s , - «^AP0 

s 



where Aw sz is the sz component of Au. We solve for pressure AP using the condition V • Au = 0: 



V • Au = V sz • u sz + V • u = V sz • u sz - kAP (^y^j = AP • k, 

At k^ = 0, Au^ = 0, this condition simply becomes V ' sz • u sz = 0. 

Therefore, at a given AT, Ap, we can solve Au, AP from this set of linear equations. 



-V S2 • u sz . 



b. Numerical computation of linear response matrix 



To solve for stability, we discretize the system in the s, z directions in the same way that we did when solving for the 
equilibrium state. Because we are only solving for linear stability and the system starts off uniform in the (j) direction, 
we don't need to discretize in the (j) direction. For perturbations at a certain wavenumber k^, the temperature 
derivative is, to first order: 
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dTi. 
dt 



. Plj . e -^=M{[T + 8T}^(p + 8p) tJ 



Metabolism 



Conduction in d> direction 



Vij f^i w 

(i> 3 >) 



• H[k ( Pij + S Pij ) , k ( Pi , jt + 5p vr )] ■ (T + 5T).,., - (T + 5T) 



Conduction in s, z directions 



v i - 



{i'j') 



Air flow in s, z directions 



Air flow in 6 direction 



Note the slight modification of the upwinding terms, where we have also included a (j) component to represent the 
influx or outflux of air in the (j) direction. 

The density derivative is, to first order: 



d PiJ c -i^ 
dt 



Jo 



^13,1 3 



(i'j') 



Pij H~ Pi' j' 



(P b + 5P b ) zi -(P b + 5P b ) irj , 



Movement in 6 direction 



Movement in s,2 directions 



The airflow is solved using the set of linear equations: 



A •/ •/ 

H (n(pij + 5 pij), K,(pi>j> + 6 p%' j')) • lJ,lJ 



(p + - (p + sp) ifjf + (z, r - z l3 ) • ( ^ + mj y + ~ 90 



SP is set such that the divergence in the <p direction negates the divergence in the s, z directions, 

J2 = 5PijK(pij)w 2 ( ^ j 

{i'j') \ tJ / 



for all cells (i, j). 



